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Abstract. Belief propagation (BP) is a message-passing method for solving probabilistic graphical models. 
It is very successful in treating disordered models (such as spin glasses) on random graphs. On the other 
hand, finite-dimensional lattice models have an abundant number of short loops, and the BP method is still 
far from being satisfactory in treating the complicated loop-induced correlations in these systems. Here we 
propose a loop-corrected BP method to take into account the effect of short loops in lattice spin models. 
We demonstrate, through an application to the square-lattice Ising model, that loop-corrected BP improves 
over the naive BP method significantly. We also implement loop-corrected BP at the coarse-grained region 
graph level to further boost its performance. 

PACS. 02.70.Rr General statistical methods - 75.10.Nr Spin-glass and other random models - 07.05.Pj 
Image processing - 05.50.+q Lattice theory and statistics (Ising, Potts, etc.) 


1 Introduction 

Belief propagation (BP) is a message-passing method for 
solving probabilistic graphical models. It was developed in 
the computer science research field [1] and, independently, 
also in the statistical physics field along with the replica- 
symmetric mean field theory jjj. For spin glass physicists 
the BP method is commonly referred to as the replica- 
symmetric cavity method. The basic physical idea behind 
BP is the Bethe-Peierls approximation 0HE], which as¬ 
sumes that if a vertex is deleted from a graph, all of its 
nearest neighboring vertices will become completely un¬ 
correlated in the remaining (cavity) graph. BP has good 
quantitative predicting power if the graph’s characteris¬ 
tic loop length is much longer than the system’s typical 
correlation length. 

The BP method is exact for models defined on a tree 
graph which contains no loops. A finite-connectivity ran¬ 
dom graph contains many loops, but the typical loop length 
increases logarithmically with the total number of vertices 
in the graph, and BP also performs excellently on suffi¬ 
ciently large random-graph systems. A lot of random com¬ 
binatorial optimization problems and random-graph spin 
glass models have been successfully solved by BP and the 
replica-symmetric mean field theory during the last two 
decades [5]. 

Finite-dimensional lattice models have an abundant 
number of short loops, which cause complicated local cor¬ 
relations in the system. The correlation length of the sys- 
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tern at sufficiently low temperatures often exceeds the 
characteristic length of short loops. At the moment, BP 
is still far from being satisfactory in treating the compli¬ 
cated loop-induced local correlations in these systems. In 
recent years the generalized belief propagation (GBP), as 
a promising way of overcoming the naive BP’s shortcom¬ 
ings, has been seriously explored (7ll8ll9 UTf)lll 1) . The GBP 
method is rooted in the cluster variational method mm 
and it abandons the Bethe-Peierls approximation. 

Here we explore a simple way of improving BP while 
still keeping the Bethe-Peierls approximation. We propose 
a loop-corrected BP method to take into account the effect 
of short loops in lattice spin models. The loop-corrected 
BP method, as a hierarchical approximation scheme, is 
conceptually straightforward to understand, and its nu¬ 
merical implementation appears to be easier than the GBP 
method. As a proof of principle, we apply loop-corrected 
BP to the square-lattice Ising model for which exact re¬ 
sults are available, and demonstrate that it indeed signifi¬ 
cantly outperforms the naive BP. We also implement loop- 
corrected BP at the coarse-grained region graph level [H] 
to further boost its performance. Our numerical results on 
the square-lattice Ising model indicate that loop-corrected 
BP might be a preferred method than GBP. 

The actual applications of loop-corrected BP to the 
Edwards-Anderson spin glass model m on the square 
lattice and especially on the three-dimensional cubic lat¬ 
tice will be carried out in a follow-up paper. As poten¬ 
tial practical applications, we suggest that loop-corrected 
BP might be useful in two-dimensional image processing 
tasks, such as image recovery p6j . 
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For reason of clarity, in the remaining part of this pa¬ 
per we describe the loop-corrected BP method using the 
square lattice as a representative example. 

Let us finish the Introduction by noting that loop cor¬ 
rection to BP has been a focusing issue in the last decade 
and various protocols have been investigated Hamm 
mmmmm)- For example, the proposal of Mooij and 
co-authors {20] (see also [23]) considers the correlations 
among the neighboring vertices of a given focal vertex by 
computing the joint distribution of all these neighboring 
vertices’ spin states. This proposal abandons the Bethe- 
Peierls approximation and it is computationally expen¬ 
sive, while its performance on square-lattice spin models 
seems to be inferior to that of the conventional cluster 
variational method [20'. Another interesting approach [21] 
is based on the self-avoiding walk tree representation of a 
loopy graph [23] and its full potential is yet to be explored. 
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Fig. 1. The square lattice G with periodic boundary condi¬ 
tions. There are L (here L = 5) vertices on each boundary line, 
and the total number of vertices is N = L x L. 


2 The lattice spin system 

Let us consider a periodic square lattice G of width L 
containing N = L x L vertices, see Fig. [l] (the numerical 
results shown in Fig. [3] and Fig. [6] correspond to L = oo). 
Each vertex m £ {1, 2,..., N} of this lattice has a spin 
state <r m £ {—1, +1} and it interacts with its four nearest 
neighboring vertices. The interaction between two vertices 
m and n is represented by an edge in the lattice and this 
edge is denoted as (m, n) in our following discussions. The 
set formed by all the nearest neighboring vertices of vertex 
m is denoted as dm, i.e., dm = {n : (m,n ) £ G}. For 
the particular example of Fig. [lj dm = {l,h,n,r} and 
dn = {m,i,o, s}. In addition, we denote by dm\n the 
set obtained by deleting vertex n from the set dm, e.g., 
dm\n = {l, h, r} and dn\m = {*, o, s}. 

We denote a microscopic spin configuration of the whole 
lattice G as a, that is, a = {a i, cr 2 ,..., cttv}- The energy 
for each of the possible microscopic configurations is 
defined as 

J i3 a l a j , (1) 

i£G (i,j)£G 


subset of the 2 N microscopic configurations |2B] that are 
mutually reachable through a chain of local spin flips. 

3 The Bethe-Peierls approximation and the 
belief-propagation equation 

We now briefly review the BP method. Within a macro¬ 
scopic equilibrium state S , the marginal probability dis¬ 
tribution q m {o'm) for the spin state of a single vertex to is 
defined as 

Y' 5 a e~P E (-^ 

’ (2) 

where 5° is the Kronecker symbol such that 6° = 1 if a = 
a and 5° = 0 if a ^ cf; /3 = 1/T is the inverse temperature; 
the superscript ' of the summation symbol means that 
the summation is over all the microscopic configurations 
a belonging to the macroscopic state S. 

Since vertex to interacts only with the vertices in dm, 
we divide the total energy E(a) into two parts: 


where h,® is the local external field on vertex i, and is 
the spin coupling constant of the edge (i, j). In the limiting 
case of Jij = +J for all the edges, this model is the ferro¬ 
magnetic Ising model [25] . In the other limiting case of the 
Edwards-Anderson spin glass model, each edge coupling 
constant Jij is set to be + J or — J with equal probability 
and independently of all the other coupling constants Ba¬ 
in the numerical calculations of this paper we choose the 
energy unit to be J, which is equivalent to setting J = 1. 

Let us denote by S a macroscopic equilibrium state 
of the system at a given temperature T. When T is suf¬ 
ficiently high the system has only a single macroscopic 
state, then S contains all the 2 N microscopic configura¬ 
tions. At certain critical temperature value T c an ergodicity- 
breaking transition may occur in the configuration space 
of the system, then the system at T < T c has two or 
even many macroscopic states, each of which containing a 


E(a) = 


-h° 


- E J - 


mn'- J n 


n£dm 


+ E\ m (a\ m ) , (3) 


where E\ m (a^m) i s total energy of the cavity lattice 
G\ m formed by deleting vertex m from the original lattice 
G (see Fig. |]): 


E\m(v\ m ) = - > ( 4 ) 

*6 G\ m 

and ajm = Wj '■ J ^ G\m}- After inserting Eq. (|3]) into 
Eq. ([2]), we obtain that 

e Bh - a E <7wfehJ II 

/ >. __ _ ^drn _ nedm _ 

E q\m(<L dm ) n ( } 

<y rn <Ldm n(zdm 
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Fig. 2. The square lattice G\ m obtained by deleting vertex m 
(and all its attached edges) from the lattice G of Fig. [l] Such 
a lattice is referred to as a cavity lattice. 


where q_ dm = {a n : n £ dm} denotes a microscopic spin 
configuration of the vertices in set dm, and q\m(<ldTn) is 
the probability distribution of <Tg m within the macroscopic 
equilibrium state S of the cavity lattice G\ m : 

Q\miSSm) 

- 

Since vertex m is absent in the cavity lattice G\ m , one 
expects that the correlations among the vertices of set 
dm are much weaker in G\ m than in the original lattice 
G. Following the idea of Bethe and Peierls [30, let us ne¬ 
glecting all the remaining correlations among the vertices 
of dm. in G\ m and approximate q\m{<Ldm) by the following 
factorized form: 

9\m(<Zdm) ~ < ?n\m( cr n) ; (7) 

n£dm 


EK e 




V'- 


w n 

,~PE\ m {g _. 


n£dm u d' n 


( 6 ) 


where 9 n \ m (<r n ) is the marginal probability distribution 
of the spin state of vertex n in the cavity lattice G\ m . 
Inserting Eq. 0 into Eq. 0) we obtain the following ap¬ 
proximate expression for <?„)"( cr m ): 


q m (cr) = 


e P h m (T J] \j2e^ J ^q n \ m (a n ) 

n^dm^Cn 


n E e^J^q n \ m (a n ) 

crm n£dm L cr n 


■ ( 8 ) 


Similar to Eq. Q, we can apply the Bethe-Peierls ap¬ 
proximation on the cavity lattice G\ m to compute the 
marginal probability distribution q n \ m (a m ) of vertex n: 


Qn\m,y) 


e 0h > n [E^ J “N A n(^)' 

i£dn\m 


n E ^ J ^q An y) 

&n i£dn\m L o'i 


• (9) 


The above equation is referred to as a belief-propagation 
equation in the literature pQ. The BP equation is a self- 
consistent equation. We can iterate Eq. @ on all the edges 


of the lattice G and, if this iteration reaches a fixed point, 
then use Eq. ([8]) to compute the mean spin value of any 
given vertex m in the lattice. 

The above-mentioned mean field theory is very suc¬ 
cessful in quantitatively predicting the properties of spin 
models on random finite-connectivity graphs m- How¬ 
ever, when applied on the square-lattice Ising model with 
no external field, it predicts a transition between the para¬ 
magnetic phase and the ferromagnetic phase at the criti¬ 
cal inverse temperature /3 ss 0.3466, which is considerably 
lower than the exact value (3 C = ln(l + x/2)/2 ss 0.4407 
[28[l29]. see Fig. [3j For the Edwards-Anderson spin glass 
model on the periodic square lattice (again with no ex¬ 
ternal field), the paramagnetic solution of the BP equa¬ 
tion ^ becomes unstable as /3 exceeds certain threshold 
value p c (L) which is a decreasing function of lattice size 
L and lim^oo f3(L) ss 0.370 T3]; BP converges to a non- 
paramagnetic fixed point at /3 slightly beyond /3 C (L), but 
it fails to converge at j3 > 0.66 (see, for example, f?0]). 
These latter results are contradicting with the strong nu¬ 
merical evidence that the two- 

dimensional Edwards-Anderson model is in the paramag¬ 
netic phase at any finite /3. 

The mean-field equations ([8]) and 0 are not accurate 
in treating lattice spin models. We now develop a loop- 
corrected belief propagation numerical scheme to better 
considering the complicated effect of short loops. 


4 Loop-corrected belief-propagation equation 


We notice that, due to the abundance of short loops, the 
naive BP equations 0 and 0 generate a spurious self¬ 
field on each vertex of the lattice. By definition the proba¬ 
bility distribution q n \ m (a n ) in Eq. (8j) is completely inde¬ 
pendent of vertex m, but if we use Eq. (|9|) then q n \ m {a n ) 
will be strongly affected by m. To explain this point by 
an example, let us consider the path m-h-i-n in Fig. [2] 
9ra\m( (J n) depends on qi\ n (<Ji), which in turn depends on 
qh\i{&h), which in turn depends on < 7 m \fc,(cr m ). Similarly, 
other short paths between vertex n and vertex m will bring 
additional dependence of 9 ra \ m (cr„) on the ‘deleted’ vertex 
to. Since all the input probability distributions to vertex 
to in Eq. ([8]) actually are affected by vertex to, the result¬ 
ing marginal probability distribution g m (cr m ) contains the 
self-field of vertex to to itself. This self-field effect is not 
real but is an artifact of the naive BP equation 0. 

We need to modify Eq. ([9]) to remove this spurious 
self-field effect. Actually, if we strictly follow the Bethe- 
Peierls approximation, the expression for the probability 
distribution g n \ m (<r„) is not Eq. (|9j) but the following: 




i(zdn\m 


n E 

<y' n i(zdn\m °i 


( 10 ) 

where 9i\{ mjTl } (cq) is the marginal probability distribution 
of vertex i’s spin state in the cavity lattice G\{ m ,n} with 
both vertex to and n being deleted (see Fig. Bj). 
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Fig. 3. The inverse temperature f} c at the ferromagnetic phase 
transition point of the square-lattice Ising model (no external 
held). The results obtained by the belief-propagation equation 
(BP, plus symbols) and those obtained by the loop-corrected 
belief-propagation equation (LC-BP) with memory capacity 
C = 2 (star symbols) and memory capacity C = 3 (cross sym¬ 
bols) are compared with the exact value /3 C fts 0.4407 (marked 
by the horizontal dashed line). Each square region of BP and 
LC-BP contains nxn vertices, with n being the number of ver¬ 
tices along one boundary line of the square region. We can fit 
the data by the function /3 C = /3%° — cn \ with /3^° = 0.4490, 
c = 0.1109 and 7 = 0.6075 (for BP, bottom dashed curve), 
0“ = 0.4421, c = 0.0746 and 7 = 0.8357 (for LC-BP at 
C = 2 , middle solid curve), and /3£° = 0.4417, c = 0.0517 
and 7 = 0.8071 (for LC-BP at C = 3, top dotted curve). 


In general, for any given vertex set <j> and a vertex n 
that is adjacent to at least one vertex in this set <f>, we 
denote by q n \^(<J n ) the marginal probability distribution 
of vertex n ’s spin state in the cavity lattice obtained 
by deleting all the vertices of <f> from the original lattice G. 
Under the Bethe-Peierls approximation, this probability 
distribution can be determined through 
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Fig. 4. The cavity square lattice obtained by deleting 

vertices m and n (and all the attached edges) from the lattice 
G of Fig. [l] 


vertex n can contain at most two vertices (i.e., memory 
capacity C = 2). Under this additional restriction, then 
for the two vertices l and r in Fig. [ 4 ] we have 


Ql\{m,n} ( 07 ) CK 

eP h ?' 

"E 

{,0 a l Jig' 

<*9^ 

9 Qg\{l 1 m}(°'g)] 



X 

E« 

gflcriJlk&k 

Qk\{l,m} (^/c)] 



X 

E a i 

g/3<70 JlqCTq 

^g\{i,m}( cr g)] j 

(12a) 

Qr\{m,n} (^"r) 

eP h °r- 

CTr E 

„0a r J rc 

1 q Qq\{r,m} (^g)] 



X 

E- 

e /3o> J rw o 

Qw\{r,m } (^u;)] 



X 

E- 

e Pa T J rs o. 

(Zs\{r,n} (^s)] • 

(12b) 


We consider g s \( r . n \(a s ) instead of <? s \{i-,m}(o' s ) in the last 
line of Eq. (12b I because vertex n has stronger influence 
to vertex s than vertex to. The probability distribution 
Qs\{r.n} (°"«) °f Eq. (12b) can be computed through 




i£dn\c/> L cr; 


£e^« n 


i£dn\4> L a 


e^" Jniai qi\^ n y(ai) 


( 11 ) 

where dn\<j> = dn — <pndn denotes the vertex set obtained 
by deleting all the vertices of dn that are also belonging 
to set <j), and {cj), n} = <f> U {n} is the vertex set obtained 
by adding vertex n to set (f >. 

Equations Q, ( |To| ) and 0 form a hierarchical se¬ 
ries of self-consistent equations and we refer them collec¬ 
tively as the loop-corrected belief-propagation equation. 
For practical applications we have to make a cutoff to this 
message-passing hierarchy, so that a closed set of equa¬ 
tions can be obtained and can be iterated numerically. 

In the remaining part of this paper we mainly consider 
the simplest nontrivial cutoff by requiring that the vertex 
set cf> of the cavity probability distribution q n \(f, of any 


9a\{„,r}K) E^'^AI-mM 

• (I 3 ) 


When we apply Eqs. (12) and (13) to the square-lattice 
Ising model, we obtain a critical inverse temperature /3 C ~ 
0.3716 for the ferromagnetic phase transition, which is 
considerably better than the prediction of the naive BP, 
see Fig. [3] This is an encouraging result. We can further 
improve the performance of the loop-corrected BP mean 
field theory by allowing the set </> of deleted vertices in 
Eq. 0 to contain three or even more vertices. For ex¬ 
ample if the memory capacity is set to C = 3 the value of 
f3 c estimated for the ferromagnetic Ising model increases 
to /3 C « 0.3896 (see Fig. [ 3 ]). 

The mean magnetization (<r m ) of vertex to. and the 
mean spin correlation (a m a n ) between vertex to. and n 
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are estimated through the following equations: 

(&m) — ^ } (Ida) 

(7 m 

E cr mO’ rl e / 3 Jm '* ,T ™ CT "g m \ n (CT m )g n \ m (<j„) 

( ' rTm<Tn ' ) E e /3Jm ’ lCTmCT ’ , 9m\n(c r m)gn\ m (c r n) 

<J m ,ff n 

(14b) 

The mean energy of the whole system is then 
N 

(#) = ~ E h 0 m {(Jm) - E Jrnn{Vm<Jn) ■ (15) 

m—1 (m,n)(£G 

(E ) of course depends on the inverse temperature /3, let 
us emphasize this dependence by (E)p. The free energy 
F(/3) of the system is related to the mean energy through 
(E)p = ^§p-, namely 



Fig. 5. The square lattice coarse-grained as a region graph TZ. 
Each square region ( 71 , 72 ,...) contains n x n vertices and all 
the interactions within these vertices (n = 2 in this particular 
example). Two nearest-neighboring regions interact through n 
edges. 


P 

F{P) = J (E)p'&fi — -^7VTn2 . (16) 

o 


5 Loop-corrected belief propagation at the 
region graph level 


In essence, the loop-corrected BP mean field theory of the 
preceding section tries to completely eliminate the effect 
of a deleted vertex m to the cavity lattice G\ m through the 
BP hierarchy Eqs. (10 1 and ©• But the loop-corrected 
BP hierarchy is alsobased on the Bethe-Peierls approx¬ 
imation and it does not consider any of the short-range 
correlations that are discarded from this approximation 
(e.g., the correlations among the vertices I, h, n , and m in 
the cavity graph G\ m of Fig. f2|) . To take into account more 
short-range correlations, we iollow the work of Zhou and 
Wang HH and construct the loop-corrected BP equation 
at the coarse-grained region graph level. 

In the example of the square lattice, we completely 
cover the vertices of the whole lattice by a set of square 
regions without any overlap between the regions. Each 
square region contains nxn vertices and all the interaction 
edges within these vertices, see Fig. [5] Two neighboring 
regions interact with each other through the n edges in 
between, and they are therefore considered as being con¬ 
nected at the region level. The region graph TZ constructed 
in this way, with each vertex representing a local square 
domain of nxn vertices, has the same topology as the 
original square lattice G. 

The loop-corrected BP hierarchy can then be obtained 
for this region graph TZ. Consider the region 75 of Fig. [5] 
as an example. Let us define q^ 5 \ 72 (a m ,a n ) as the prob¬ 
ability of vertex m taking spin value a m and vertex n 
taking spin value a n in the cavity region graph TZ\ l2 ob¬ 
tained by deleting region 72 from TZ. Other joint proba¬ 
bility distributions can be defined in a similar way, e.g., 


< 7 75 \{ 7 lj 72 }(cr m , er„) is the joint probability distribution of 
<r m and a n in the cavity region graph 7£\{ 71l72 } (with re¬ 
gions 71 and 72 being deleted). If we restrict the set </> of 
deleted regions in memory to containing two regions at 
most (i.e., memory capacity C = 2), we obtain that 

^75X72 (pm, &n) 'y " ^ 75 

x \52*„a 9 e- fiK '™ 9 74 \f 72 , 7B } (<7z, Tq )] 

X ^ 7875 9 78 \{ 7 2,7 5 }( <T '!«: a x)\ 

x 7675 976 X 172 , 75 }(®t> o’o)] ) (17a) 

975 X 171 , 72 ! { a mi O’n) 7^ ay,a- a e ^ 75 

X E *, ,<r q 3 74 \ {71 , 75 } , <J q )\ 

X ^ 7875 9 7 8\{ 7 2,75} (.^w: *Tc)] 

x E— e -^-^\ {72 , 75 >(— 0 )] > ( 17b ) 

975 X 172 , 74 } { a mi O’n) CX ^~^ q- r ,q- a e ^ T5 

X 7 ' iy w ,r7x e ^ 7875 9 7s \1 74 ,75} {p'wt (7 x)] 

X E — 0 e_/3£7675 9 7 e\l72,75}(^Wo)] ■ (17c) 

In the above expressions, the quantity E 1 denotes the in¬ 
ternal energy of a region 7 , for example 

-^75 &s) — h r (T r h s (J s 

Jmn®m®n Jns®n®s Jrs&r&s Jmr^m^r •> (18) 

and E 77 / is the interaction energy between region 7 and 
region 7 ', for example 

-£7475 (^7 5 &rm &r) — Jqr&q&r • ( 19 ) 

As Eq. p~7| demonstrates, all the correlations within 
each region are precisely considered by summing over all 
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the 2 ra microscopic configurations of this region. In the 
practical implementation, the internal state summation is 
achieved through a numerical scheme that is efficient both 
in terms of computing time and in terms of needed mem¬ 
ory (see Appendix A for details). By increasing the region 
size n we can include more and more short-range correla¬ 
tions and achieve more precise quantitative predictions. 

For the two-dimensional Ising model we have compared 
in Fig. [3] the results obtained by the conventional region- 
graph BP of [21 and those obtained by the present region- 
graph loop-corrected BP. When the memory capacity is 
set to C = 2 (the smallest nontrivial value), the iteration 
process of loop-corrected BP demands the same order of 
computational cost as that of BP, yet at each value of the 
square-region size n the improvement of loop-corrected 
BP over BP is always significant, suggesting that loop- 
corrected BP is a much better choice than the naive BP for 
treating finite-dimensional lattice systems. Figure [6] com¬ 
pares the exact spontaneous magnetization of the square- 
lattice Ising model with the predictions obtained by BP 
and LC-BP (<7 = 2). At each value of the region sizes used 
(n = 1, n = 3, or n = 5) the improvement of LC-BP over 
BP is again significant. 

It also appears that loop-corrected BP (with memory 
capacity C = 2) outperforms the GBP method of Yedidia 
and coworkers [7|. When the square-region size is set to 
n = 2, GBP predicts the critical inverse temperature of 
ferromagnetic phase transition to be /3 C « 0.4126 [TO]; a 
slightly better result is achieved by the loop-corrected BP 
method at square-region size n! = 2n = 4, which reports 
a value of /3 C ~ 0.4186. The GBP with square-region size 
n = 4 predicts a value of /3 C ss 0.429 EH; this result is 
marched by the loop-corrected BP at square-region size 
n' = 2 n = 8, which reports a value of /3 C ~ 0.4290. We 
might therefore conjecture that GBP at square-region size 
n and loop-corrected BP at square-region size n! = 2 n 
have comparable prediction power. Under such an assump¬ 
tion we can then argue that loop-corrected BP will be a 
better choice than GBP: (1) the iteration process of GBP 
is much more complicated than that of loop-corrected BP; 
and (2) the required computer storage space of a GBP 
message is of order 0(2 n / 2 ), making it unpractical to 
set the square-region size n > 6; (3) the required stor¬ 
age space of a loop-corrected BP message is only of order 
0(2 n ), so we can set the square-region size to n = 20 or 
even larger values. It should be pointed out that good per¬ 
formance of GBP can be achieved by increasing the size 
of the largest region one-dimensionally rather than two- 
dimensionally (see [39] and 0). It will be helpful to per¬ 
form a comparative study by implementing LC-BP also in 
such a non-symmetric way. We leave this point for future 
investigations. 

We can further improve the performance of the loop- 
corrected BP method by increasing the memory capacity 
C (but at the cost of introducing many more cavity mes¬ 
sages, see Appendix B). For the square-lattice Ising model, 
the results obtained by loop-corrected BP at C = 3 are 
also shown in Fig. [3] to compare with the results obtained 
at C = 2. We find that increasing C from C = 2 to C = 3 



P 


Fig. 6. The spontaneous magnetization (the mean spin value 
(a) of a vertex) of the square-lattice Ising model. The results 
obtained by belief propagation (BP) and those obtained by 
loop-corrected belief propagation (LC-BP) with memory ca¬ 
pacity C = 2 are shown together with the exact results (the 
solid line). Each square region of BP and LC-BP contains nxn 
vertices, with n = 1, 3, or 5. 

does not bring a dramatic improvement to the prediction 
of /3 C . Considering the high computation cost required for 
C > 3 (see Appendix B), if higher numerical precision is 
needed, it is more practical to increase the square-region 
size n but to keep the memory capacity at (7 = 2. 


6 Conclusion 

To summarize briefly, in this paper we described the main 
ideas of the loop-corrected belief propagation method and 
carried out an initial performance test on the square-lattice 
Ising model. The results in Fig.[3]and Fig.[6]clearly demon¬ 
strate that loop-corrected BP with memory capacity C = 
2 is much superior to the naive BP method, which is equiv¬ 
alent to loop-corrected BP with memory capacity (7=1. 
The performance of loop-corrected BP further improves as 
the memory capacity is increased to C = 3 or even larger 
values. 

Our numerical results on the square-lattice Ising model 
also indicate that, compared to the generalized belief prop¬ 
agation method of Yedidia et al. [7j, the loop-corrected BP 
method (simply with memory capacity (7 = 2) can achieve 
the same or even higher level of precision at much re¬ 
duced computation cost. In addition, we wish to point out 
another very important advantage of the loop-corrected 
BP method: just as the survey propagation method is a 
natural extension of the naive BP method [61127] . follow¬ 
ing the discussion of El we might extend loop-corrected 
BP into the loop-corrected survey propagation method to 
study disordered lattice models in the low-temperature 
spin glass phase, where ergodicity of the configuration 
space is broken. 












Hai-Jun Zhou and Wei-Mou Zheng: Loop-corrected belief propagation for lattice spin models 


7 


For the loop-corrected BP method really to be a help¬ 
ful tool, it should be capable of giving good quantitative 
predictions on single instances of disordered lattice mod¬ 
els. The performance of loop-corrected BP on the square- 
lattice and cubic-lattice spin glass models will be investi¬ 
gated and be reported in a forthcoming paper. 
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Appendix A: Message updating for a square 
region 


To perform region-graph BP or loop-corrected BP itera¬ 
tion on a square lattce, the most demanding task is com¬ 
puting the joint probability distribution of spin states for 
the vertices on the boundary of a region. Let us con¬ 
sider the concrete example shown in Fig. [TJ The central 
(C) square region contains n x n vertices with n = 6 , 
and it receives messages from three other square regions 
on the left (L), bottom (B), and right (R) side. Denote 
by g_ T = ((T 2 , 03 ,..., ay) a generic spin configuration for 
the n vertices on the top (T) boundary of the central re¬ 
gion. This spin configuration is affected by the interactions 
within the central region and the interactions between the 
central region and the three neighboring regions, and its 
probability distribution Pt{q_t) is expressed as 


Pt{<L t ) oc e P Ec ^2P L (a L )e P El ’ c 

—C\T °L 

X [j2 p B^B)e~ 0EB - c ] [J2 p R^r)^ Er ’ c 

—B 2.R 


( 20 ) 


In this expression, a c \ T = ( 027 , oyg) • • •, 055 , 056 ) is a spin 
configuration for all the other (n — 1) x n vertices of the 
central region except the n vertices at the top boundary, 
and Ec is the total internal energy of this central region; 
g_ L = (oq, 026 , ■ • •, 022 ) is a spin configuration for the n 
boundary vertices of the left region, and Pl(<Ll) is an 
input probability distribution of q_ L , and E^c is the in- 
teration energy between the left and the central region; 
similarly, <r R = ( 0 - 33 , 012 , • 0 ’s) is a spin configuration 

for the boundary vertices of the right region, Pr{q_ r ) is 
an input probability distribution of a R , E Rt c is the in¬ 
teraction energy between the right and the central region, 
and a B = ( 020 , 019 , ■ ■ ■, 0 15 ) is a spin configuration for the 
boundary vertices of the bottom region, P B (a B ) is an in¬ 
put probability distribution of a B , E B ^c is the interaction 
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Fig. 7. The central square region contains n x n vertices 
(n = 6) and it interacts with the three neighboring square 
regions (partly shown) on the left, bottom, and right side. 


energy between the bottom and the central region. Notice 
that the LC-BP equations (12), (13), and ( fl7| ) all have the 
same form of Eq. (20 1 . 


According to Eq. (20 1 , one needs to sum over a total 
number of 2 "(” +2 ) different spin configurations to deter¬ 
mine the output probability Pt{qLt) of a single spin con¬ 
figuration a T . A naive application of Eq. (20) is therefore 
feasible only for very small values of n (e.g., n < 3). 


We now introduce a numerical trick that greatly ac¬ 
celerate this summation process. By this simple trick we 
reduce the total number of needed operations to sum over 
all the spin configurations from 0 ( 2 n ( n + 2 )) to 0(n 2 2 n ), 
and also dramatically reduce the total amount of storage 
space needed in the numerical computation. 


First we notice that, due to the binary nature of the 
spins, a generic probability distribution p(oi, oy, • • ■, <J n ) 
over n spins can be written in the following form: 


11 1 

P(oi, • • •, oh) = c ^...s n ^ S i^2 

si—0 S 2 —0 s n —0 

( 21 ) 

where s, £ {0,1} for i = 1, 2,..., n and {c SlS2 ... s „} is a set 
of 2 ra coefficients, with C00...0 = 2 _n due to the normal¬ 
ization constraint. Therefore the probability distribution 
p(a 1 , ay,... , cr n ) is completely characterized by the coef¬ 
ficient set {c SlS2 ... s J. 

Due to the fact that 


e 0 h i<Ti s cosh (/3hi) [l + tanh(/3hi)ai] , (22a) 

e P J a' cria i' = cosh(/3Jij/) [l + tanh(/3Jjj/)a-j0-j/] , (22b) 


















































































Hai-Jun Zhou and Wei-Mou Zheng: Loop-corrected belief propagation for lattice spin models 


then for i, j e {1,2,..., n} (i < j) and i' ^ {1,2,..., n}, 



...,a n ) = cosh(/3 hi) ^ 0{V s 2 2 ■ 

S\S 2 ...Sn 


[ C SlS2- 

..s n , 

(23a) 


ai<Ti 'p((T U ...,(j n ) = 2 cosh (fiJu') 


X 

M 

-S 1 - s i~l -Si Si+l S n 

G l • • ' a i-l G i' a i+l • • • s n 


s 1 s 2 — s 

n 


x [(1 

H - Cs 1 s 2 ...s n 5 

(23b) 


ePJijOiajp^ i ) ..., a n ) = cosh(/3Jy) x 

E ^l 1 ^ 2 ■ ■ • S n l [ C 'S 1 s 2 ...s t , - I - 

s 1 s 2 ...s n 


ta,nh(/3Jij)c sl ... Si _ lSiSi+1 ... Sj _ lS:i g 3 . +1 ... s „] , (23c) 
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Fig. 8 . Given an input probability distribution 
Pb{ct 2 o, ■ ■ ■ , C 15 ) for the set {20,19,, 15} of vertices on the 
bottom row, the probability distribution Qb(<74i, • • •, 036 ) 
for the set {41,53,55,56,50,36} of boundary vertices can be 
determined by recursion from the bottom row up to the top 
row. 


{2,28,46,47,31,7} through the following expression: 


where s,; = 1 if s,; = 0 and Si = 0 if Sj = 1. Equation 
(23) therefore gives a set of rules on how the coefficients 
set {c SlS2 ... Sri } changes as p(<7i,..., cr n ) is perturbed by 
multiplication and summation. 


We simplify the computation of Eq. (20) by treating 
the three input probability distributions separately. For 
example, starting from the input probability distribution 
Pb(&2 0, 019, ■ ■ ■, 015) of the bottom region (see Fig. [ 8 ]), 
we obtain a probability distribution Qb(& 41 , <753, ..., 1736) 
for the set of n boundary vertices {41, 53, 55, 56, 50, 36} 
through the following recursive process: ( 1 ) initialize the 
coefficients set of Qb(-) to be identical to that of Pb(-); (2) 
then consider sequentially all the n vertical edges (20,41), 
(19,40), ..., (15,36) between the central and the bottom 
region and m odify the coefficients set of Qb(-) according 
to Eq. (23b); (3) then consider sequentially all the (n— 1) 
horizontal edges (41,40), (40,39), ..., (37,36) between the 
set of vertices {41,40, 39, 38, 37, 36} and further modify 
the coefficients set of Qb(') according to Eq. (23c); (4) 
then consider sequentially all the (n — 1 ) external fields 
on the set of internal vertices {40, 39,..., 37} and f urthe r 
modify the coefficients set of Qb(‘) according to Eq. (23a); 
(5) repeat the previous three steps on the row co ntain- 
ing the set of vertices {53,52,51,50}: apply Eq. ( 23b| 
on the set of ve rtica l edges {(40, 53),..., (37, 50)} am 
then apply Eq. (23c I on the horizontal ed ges (53,52), 
(52,51) and (51, 50)}, and then apply Eq. (|2.3a[ ) on the 
internal vertices 52 and 51; ( 6 ) finally, apply Eq. ( |23b | 
on the edges (52,55) and (51,56), and apply Eq. ( |23c | 
on edge (55,56) and then output the coefficients set of 
Qb(c 41 , 053 , 055 , 056 , 050 , 036 )- 

The joint probability distributions Ql(p 2 ,..., (741) for 
the set of vertices {2, 28,46,55,53,41} and Qr((j 36,..., ( 77 ) 
for the set of vertices {36, 50, 56,47, 31, 7}, see Fig. [7j are 
obtained through the same recursive process starting from 
Pl(-) and Pr(-), respectively. The only additional feature 
is that we now need to consider the external fields of all 
the vertices i n the se two boundary sets (again through 
applying Eq. (23a) to Ql{~) and Qr(-) repeatedly). 

With these preparations, we then obtain a joint prob¬ 
ability distribution Q(ct 2 , ■ ■ ■, 07 ) for the set of vertices 


< 9 ( 02 , 028 , 046 , 047 , 031 , 07 ) OC 47^46^47 x 

E E <9l(°’ 2,028, 046, 055, 053, 04l) 

<741 ,<753><755 <756 ,<750 ><736 

X<9b(041,053, 055, 056, 050, 036) 

X <9i?(036i 050 , 0.56 , 047 , 031, 07 ) 


p/3 >^46,47 <746 <747 


E 


rr s 2 fjS 28 ^.546 a s A7 S 31 s 7 

°2 a 28 °46 °47 °31 °7 


S 2 S 28 S 46 S 47 S 31 S 7 
(L) 


E E 

S55S53S41 S 36S50«56 

„(B) 


S2 528^46^55S53S41 


~(fl) 


(24) 


In the above expression, the coefficient sets {c^/.. S41 }, 

{cif/...^}, and {cifl.. S7 } correspond to Ql(-), <9b(-), and 
Qr(-), respectively. The effect of the multiplication term 
e ,8,746,47<746<747 the coefficient set of the probabili ty d is- 
tribution Q(-) can again be obtained through Eq. (23c). 

Finally, the probability distribution Pt( 02 , • • ■, 07 ) for 
the set {2,3,4, 5,6, 7} of vertices at the top boundary is 
determined from < 9 ( 02 , 028 , 046 , 047 , 031 , 07 ) through the 
following recursive process (see Fig. [ 9 ]): (1) set the coeffi¬ 
cients set of Pt(-) to be identical to that of <9(-); (2) then 
consider the vertical edges (29,46) and (30,47) sequen¬ 
tially a nd m odify the coefficients set of Pt(-) according 
to Eq. (23b); (3) then consider all the horizontal edges 
(28, 29), (29, 30), and (30,31) between the set of vertices 
{28,29,30,31} and the external fields on vertices 29 and 
30 and fur ther modify th e coe fficients set of Pt{-) accord¬ 
ing to Eq. (|23c|) and Eq. (23a I, respectively; (4) repeat the 


operations of steps (2) and (3) on the vertical edges be¬ 
tween the top and the second row of Fig. [9j the horizontal 
edges of the top row, and the set of vertices {3,4,5,6}. 
We then output the resulting coefficient set of Pt(-) as 
the result of original computing task Eq. (201. 


It is straightforward to extend the numerical trick of 
this appendix to other values of even n and also to the 
cases of n being odd. For studying lattice models on a 
three-dimensional cubic lattice, this same trick can be ap¬ 
plied to a cubic region containing n x n x n vertices. 
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Fig. 9. Given an input joint probability distribution 
Q(<J 2 , ..., ay) for the set {2,28,46,47,31,7} of vertices on 
the bottom boundary, the joint probability distribution 
Pt{( 72 , ■ ■ ■, 07 ) for the set {2, 3,4, 5, 6, 7} on the top row can 
be determined recursively from the bottom row up to the top 
row. 
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Fig. 10. When the memory capacity is set to C = 3, each 
focal vertex/region (denoted by a filled small square) needs 
to remember the positions of the other three deleted ver¬ 
tices/regions (denoted by three unfilled small squares). In total 
we need to distinguish 29 different patterns of the three deleted 
vertices/regions, which are indexed as 00, Ola and 016, ..., 14, 
and 15. The small arrows indicate the cavity message of the 
focal vertex/region to the deleted vertices/regions. For the pur¬ 
pose of clarity we separate different patterns through the thin 
dashed lines. 


Appendix B: Loop-corrected belief propaga¬ 
tion with memory capacity C — 3 


When the memory capacity is set to C = 3, then with 
respective to a focal vertex or region (denoted by a filled 
small square in each block of Fig. 10), we need to consider 
29 different patterns of the three deleted vertices or regions 
(denoted by three unfilled small squares in each block of 
Fig. [lO]). These 29 patterns are indexed as 00, 01a and 016, 
02a and 026, ..., 13a and 136, 14, and 15 in Fig. [To| for 
the convenience of discussion. The patterns 01a and 016 
(and similarly 02a and 026, 03a and 036, ...) are related 
by a mirror symmetry. 

Each pattern of Fig. [10] is associated with a cavity 
message. For example, suppose vertices l, m, n are deleted 
from the graph G of Fig. [T] then the pattern 01a of Fig. [To] 
corresponds to the cavity message q s \{i } m,n}( IJ s) from ver¬ 
tex s to vertex n, while pattern 0 l 6 corresponds to the 
cavity message q q \u t m,n} ( 07 ) from vertex q to vertex l. 



Fig. 11. Diagram showing how the cavity message of all the 
29 patterns in Fig. 10 are iteratively determined (see main text 
for more details). For reason of clarity, for each pair of mirror 
patterns (say 01a and 016) we only draw the input edges to 
one of the patterns (01a) but not to the mirror pattern (016). 
The edges to each mirror pattern can be easily constructed 
by symmetry considerations. For example, since pattern 01a 
receives edges from patterns 026, 03a and 04a, then pattern 
016 must receive edges from patterns 02a, 036 and 046. 


As another example at the region graph level, suppose re¬ 
gions 72 , 76 and 79 are deleted from the region graph 1Z of 
Fig- E then the pattern 036 of Fig. pH] corresponds to the 
cavity message g 75 \{ 72 , 76 , 79 }(o'm, a n ) from region 75 to 72 . 

The iteration of the 29 cavity messages for the 29 pat¬ 
terns of Fig. 10 is carried out following the updating dia¬ 
gram of Fig. Each directed edge pi —>• P 2 in this dia¬ 
gram points from one pattern (say p± = 04a) to another 
pattern (say p -2 = 01 a), and it means that the cavity mes¬ 
sage of pattern P 2 is determined (partly) from the cavity 
message of pattern p\. For example, there are three di¬ 
rected edges (from patterns 01 a, 016 and 00 , respectively) 
to pattern 00 , meaning that the output cavity message of 
pattern 00 can be computed based on three inputing cav¬ 
ity messages from patterns 00, 01a and 016. In the specific 
case of Fig. |T] we have 

qr\{c,h,m}(o7r} 05 ^ r r ) aQ 6^ qr q r q q \{h,m,r} (^q)] 

X |~ ^ ' cr w 6-^ wr w r qw\{h,m,r} (^u;)] 


E< 


p/5 Jsr<7 s&t 


qs\{h,m,r} (^s)] • 


(25) 


The updating equations for the other 28 cavity messages 
can be written down in a similar way according to Fig. El 
Notice that in Fig. 11 we only draw the input edges to pat¬ 
terns 00, 14, 15 and patterns 01a, 02a, ..., 13a but not the 
input edges to all the mirror patterns 016, 026, ..., 136 to 
avoid the diagram being too complicated. We can easily 
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construct all the missing directed edges by symmetry con¬ 
siderations. For example, since pattern 03a receives edges 
from patterns 07a and 086, then pattern 036 must receive 
edges from patterns 076 and 08a. In the specific case of 
regions 72 , 76 , and 79 being deleted from Fig. [5j we have 

<?75\{72,7e,7 9 }(<r m ,*n) « 

* E tE ^+^^^ 97A{72 ,7,7,}^,^)] 

&l iOq 

X [E g 78 \ { 72 i 75 i 79 } K,^)] . 

ow : crx 

(26) 
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